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Abstract 

An  investigation  is  conducted  into  low  orbiting  satellites  about  the 
planet  Venus  with  drag  limited  lifetimes.  It  is  possible  to  specify  combi¬ 
nations  of  orbital  elements  which  result  in  orbital  lifetimes  greater  than 
some  desired  value.  These  combinations  can  be  assembled  into  Critical 
Curves  and  Critical  Surfaces.  Critical  Curves  and  Critical  Surfaces  are 
defined  as  the  curve  or  surface  in  orbital  element  space  above  which 
initial  element  sets  will  result  in  orbits  that  meet  or  exceed  lifetime 
requirements 

A  numerical  method  is  implemented  for  finding  these  combinations, 
and  three  Critical  Surfaces  are  examined.  A  "decay  threshold"  is 
selected  for  bounding  satellite  lifetime.  Numerical  simulations  of  orbital 
behavior  are  conducted,  and  polynomials  describing  the  Critical  Curves 
are  produced  for  five  different  eccentricities  at  each  of  three  decay 
thresholds.  Comparisons  of  the  effects  of  the  different  perturbations 
considered  (geopotential,  atmospheric  drag,  third  body  effects,  and  solar 
radiation  pressure)  and  of  decay  threshold  altitude  variations  are  made. 
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MINIMUM  ORBITS  ABOUT  THE  PLANET  VENUS 


I.  Introduction 


Background 

Orbital  elements  are  initially  determined  to  meet  specific  mission 
objectives.  These  objectives  may  be  in  terms  of  ground  track  coverage, 
communications  coverage,  or  other  criteria.  For  example,  Cain  [3]  has 
determined  element  sets  which  will  result  in  orbits  containing  arcs  of 
minimum  altitude  variation.  These  orbits  would  be  useful  for  missions 
utilizing  sensors  which  are  extremely  sensitive  to  altitude  variations,  hu: 
are  also  extremely  susceptible  to  decay  due  to  drag,  because  of  low 
altitudes  and  eccentricities.  All  orbital  elements  will  tend  to  stray  from 
their  nominal  values.  These  perturbations  are  caused  by  various  forces 
acting  on  the  satellite,  and  can  result  in  orbital  decay  to  the  extent  that 
the  satellite  reenters  the  atmosphere.  Satellite  lifetime  is  therefore  an 
important  concern  to  mission  planners. 

It  is  possible  to  determine  combinations  of  orbital  elements  which 
specify  ar.  orbit  which  meets  a  given  criteria  for  orbit  lifetime.  This 
could  be  described  as  a  minimum  or  critical  orbit.  If  these  critical 
element  sets  are  known  over  a  wide  range  of  element  values,  a  "Critical 
Surface"  may  be  constructed  (in  six  dimensional  element  space),  above 
which  satellite  lifetime  can  be  guaranteed  to  meet  or  exceed  the  lifetime 
criteria.  This  Critical  Surface  may  then  be  used  as  an  aid  in  determin¬ 
ing  initial  element  sets  which  would  fulfill  other  requirements  and  meet 
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a  minimum  lifetime  specification.  This  would  be  particularly  useful  for 
low  altitude  lightsats  with  limited  or  no  station -keeping  capabilities. 

Drag  is  the  predominant  perturbation  involved  in  the  decay  of  low 
altitude  satellites.  In  determining  these  minimum  orbits,  drag  therefore 
must  be  taken  into  account.  To  accurately  model  the  satellite  lifetime, 
other  significant  perturbations  must  also  be  included.  A  decay  thresh¬ 
old  must  be  selected  below  which  the  satellite  is  considered  to  have 
"reentered"  in  order  to  definitively  bound  lifetime.  In  order  to  produce 
the  Critical  Surface,  a  wide  range  of  initial  elements  is  required,  but 
oititi.  where  drag  is  not  the  most  influential  perturbation  are  not  con¬ 
sidered.  In  these  cases,  the  orbital  elements  will  be  perturbed  from  the 
nominal  valuer,  but  the  perturbations  will  not  necessarily  result  in 
cr  Ida',  decay.  Drag  is  most  predominant  for  orbits  with  low  eccentrici¬ 
ties  sc  only  orbits  with  fairly  low  eccentricities  will  be  addressed. 

Objective 

A  method  of  determining  the  Critical  Surface  was  derived  (see 
Chapter  III,  Search  Method).  Three  Critical  Surfaces  for  the  planet 
Venus,  based  on  three  slightly  different  decay  thresholds,  were  speci¬ 
fied  by  this  method.  Geopotential  effects,  atmospheric  drag,  solar  third 
body  effects,  and  solar  radiation  pressure  were  considered. 

Approach 

A  search  method  was  employed  to  determine  the  minimum  semi-major 
axis  for  a  given  eccentricity  and  inclination  such  that  the  orbit  did  not 
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decay  to  the  point  of  reentry  before  90  earth  days  had  elapsed.  If  an 
axisymmetric  planet  is  considered,  the  longitude  of  ascending  node, 
argument  of  periapsis  and  true  anomaly  may  be  set  to  any  arbitrary 
value  without  changing  the  effect  of  the  perturbations  considered,  and 

i 

the  convenient  value  of  zero  was  used.  This  limited  the  effective  orbital 
element  space  to  three  dimensions.  The  method  was  employed  over  a 
range  of  eccentricities  and  inclinations.  Polynomials  were  then  fitted  to 

» 

these  initial  element  sets  to  obtain  equations  describing  general  relation¬ 
ships  between  the  elements  in  question.  Three  Critical  Surfaces  were 
generated  in  this  fashion. 

> 

> 

> 
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n.  Analytical  Model 


Satellite  Motion 

A  satellite  affected  only  by  the  gravitational  attraction  of  a  spheri¬ 
cal,  homogeneous  planet  would  exhibit  ’’two-body,"  or  Keplerian,  motion. 
The  two-body  solution  positions  the  satellite  by  means  of  six  variables 
(semi-major  axis  a,  eccentricity  e,  inclination  i,  longitude  of  ascending 
node  n,  argument  of  periapsis  u>,  and  true  anomaly  v).  With  the  excep¬ 
tion  of  the  true  anomaly,  which  locates  the  satellite  along  the  orbit, 
these  would  not  change  with  time.  The  satellite  would  continue  to  follow 
the  same  orbit  indefinitely.  In  reality,  a  host  of  lesser  forces  act  on 
the  satellite,  "perturbing"  the  orbit  and  causing  the  orbital  elements  to 
change  over  time.  In  order  to  study  the  effects  of  the  perturbations  it 
is  necessary  to  model  both  the  ideal  two-body  motion  and  the  perturba¬ 
tions  cf  interest. 

Perturbations  may  give  rise  to  periodic  effects,  secular  effects,  or 
both.  Periodic  effects  result  in  no  net  change  in  the  orbital  elements, 
but  vary  the  elements  around  some  nominal  value.  Secular  effects  do 
change  the  elements  with  time,  and  are  of  primary  interest  here.  The 
effect  a  given  perturbation  may  have  on  a  specific  orbital  element  is  not 
necessarily  in  proportion  to  its  magnitude. 

The  position  and  velocity  of  a  satellite  may  be  found  from: 

f-a  (1) 
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where  f  is  the  inertial  acceleration  of  the  satellite  and  a  is  the  sum  of 
the  accelerations  caused  by  perturbative  forces  acting  on  the  satellite. 
This  chapter  will  present  the  relationships  and  equations  used  to  model 
two-body  motion  and  perturbative  accelerations  due  to  the  geopotential, 
air  drag,  third  body  effects  and  solar  radiation  pressure.  These  accel¬ 
erations  produce  perturbations  which  can  be  averaged  over  the  orbital 
period  to  give  the  net  change  in  each  element  over  one  period,  and 
these  averaged  equations  are  also  presented. 
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Two-Body  Motion.  The  two-body  motion  contribution  to  a  in 
Equation  (1)  is  the  acceleration  due  to  the  gravity  of  the  primary, 
modeled  as  a  point  mass.  Gravity  obeys  an  "inverse  square"  law,  with 
gravitational  attraction  falling  off  as  the  inverse  of  the  distance 
squared.  This  can  be  written  as: 

f  =  -ni5  (2) 

r3 

where  r  is  the  position  of  the  satellite  and  11  is  the  product  of  the  mass 
of  the  planet  (M)  and  the  universal  gravitational  constant  (G).  This  can 
be  shown  to  result  in  motion  according  to  the  equation  describing  a 
conic,  or  Kepler’s  Equation: 

r  -  (3) 

1  +  ecos(v) 

where  r  is  the  distance  from  the  center  of  the  planet  to  the  satellite  [2: 

20], 

Geopotential.  The  potential  function  that  results  in  the  right  side 
of  equation  (2)  is  V  -•  -p/r.  Taking  the  gradient  of  this  scalar  function 
results  in  the  vector  acceleration  a  referenced  in  Equation  (1): 

a  (4) 

Newton  showed  that  this  potential  can  be  extended  to  a  planet  that 
could  be  modeled  as  a  series  of  concentric  homogeneous  spherical  shells. 
Unfortunately,  real  planets  deviate  greatly  from  this  model,  severely  lim¬ 
iting  its  accuracy. 

At  this  point  it  is  necessary  to  develop  a  geopotential  function  that 


will  allow  a  non -spherical,  non-homogeneous  planet.  To  do  this,  Equa¬ 
tion  (4)  is  used  considering  V  as  an  unknown.  With  a  given  by  the 
acceleration  due  to  a  point  mass,  as  in  Equation  (2),  Equation  (4)  is 
integrated  over  the  surface  of  a  unit  sphere.  This  results  in  Poisson's 
Equation  for  gravitational  fields: 

VV  - -4nGii>  (5) 

from  which  any  V  may  be  found  as  a  function  of  position,  given 
u(.-v.y.z)  --  the  density  distribution  of  the  central  body.  Considering 
only  those  orbits  outside  the  planetary  surface  (where  .  y ,  z)  -  o) 
reduces  Equation  (5)  to  Laplace's  Equation: 

V2 1/-0  (6) 

Since  the  bodies  of  interest  are  nearly  spherical,  Equation  (6)  is 
expressed  in  spherical  coordinates: 


i  a  f  2dV\  i  a  [  dv\  i  a2 v  n 
r2ari/  a;  J  r2sme<3©lSin  JO )  r2sin20a(t>z 


(7) 


where  r  is  the  radius  to  the  point  of  interest,  ©  is  the  co-latitude,  and  0 
is  the  longitude.  This  is  a  linear  partial  differential  equation,  and  can 
be  separated  into  three  functions  of  r,  6,  and  <t>  respectively.  The  sepa¬ 
rated  function  of  0  describes  a  simple  harmonic  oscillator.  Considered 
with  the  appropriate  boundary  conditions,  its  solution  becomes: 

F(4>)  -  Ck  cos*4>  +  S*  sin  *0  (8) 


which  is  recognizable  as  Legendre's  Equation,  and  is  solv' d  by  the 
Associated  Legendre  Polynomials: 

F(9)-P‘(  cos9)  (10) 


Lastly,  the  function  of  r  yields  two  solutions,  only  one  of  which  is 
physically  meaningful.  This  leaves: 


F(r)-r") 


(H) 


where  /?.  is  the  planetary  equatorial  radius.  Equations  (8),  (10)  and 
(11)  are  combined  to  produce  an  infinite  series  solution  to  Equation  (7): 

VO  ,e.  X  X  IT  Pk„(,cosQ)(C  nkcosk$+  Sntsm*^i)  (12) 

r  n-0  k-C  \  *  •  J 


C ...  and  Snk  are  now  constants  defining  the  shape  of  the  gravita¬ 
tional  field  [7:  7],  The  lead-off  term  in  this  series  will  be  the  point 
mass  potential  term.  This  function,  with  C„*  and  Sn»  supplied  by  an 
appropriate  gravity  model  for  the  planet  in  question,  defines  the  gravi¬ 
tational  field  the  satellite  operates  in  and  may  be  employed  to  determine 
the  acceleration  on  the  satellite  due  to  the  geopotential: 

a,-7K  (13) 

Gravity  models  for  most  of  the  explored  planets  have  been  compiled  from 
empirical  data  [15]. 

If  the  primary  of  interest  can  be  modeled  as  an  axisymmetric  planet 
(a  planet  of  the  form  of  a  solid  of  rotation  about  the  polar  axis,  where 
the  ellipticity  describes  the  "out  of  roundness"  of  a  meridian,  but  paral¬ 
lels  of  latitude  are  perfect  circles)  then  Equation  (12)  can  be  simplified. 
In  this  case,  there  is  no  dependence  on  the  variable  0  and  the  index  k 
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is  fixed  at  zero.  The  S„*sin  term  goes  to  zero,  and  the  Cn*cosfc$  term 
becomes  CnCf  The  Associated  Legendre  Polynomial  P°  is  equivalent  to  the 
Legendre  Polynomial  Pn.  With  the  replacement  of  Cn0  with  -J„  Equation 
(12)  reduces  to: 

l'-(r)[,-|/-(T‘)"P-(“Se)]  (,4) 

which  is  the  geopotential  expansion  for  an  axisymmetric  planet.  The 
terms  resulting  from  this  series  are  called  "zonal  harmonics."  This 
model  will  be  concerned  with  this  expansion  through  the  sixth  order 
zonal  harmonics.  The  coefficients  (J2,J3,  etc.)  used  in  computing  the 
zonal  harmonics  for  Venus  are  given  in  Appendix  A,  Table  A-l. 

Averaged  Equations  for  Geopotential.  Merson  [11]  produced  aver¬ 
aged  equations  for  the  change  in  each  of  the  elements  due  to  geopoten¬ 
tial  expansion  terms  for  J z- J t  over  one  orbital  period,  and  his  results 
are  presented  here.  If  the  inclination  or  the  eccentricity  is  zero,  the 
equation  for  the  change  in  argument  of  periapsis  is  not  valid,  and  if  the 
inclination  is  zero  the  equation  for  longitude  of  ascending  node  is  not 
valid.  Merson  also  included  terms  for  the  change  in  the  elements  due  to 
( ? ; ,  since  for  earth  (J2)z  is  of  the  same  order  of  magnitude  as  J3-Jt 
This  is  not  the  case  with  Venus,  where  the  J 2  term  is  not  nearly  as 
predominant  as  with  the  earth  (see  Table  1).  Therefore,  the  (J2)z  terms 
have  been  culled  from  the  equations  presented  here. 

The  change  in  the  semi-major  axis  is  entirely  dependent  on  the 
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(J 2)2  term  and  is  therefore  negligible.  The  change  in  the  other  ele¬ 
ments  is  given  by: 


Ax  =  2ji 


/?. 

a(l  -e2) 


(15) 


where  Ax  is  replaced  with  Ae,  Ai,  Afl  or  Au>,  and  xn  is  replaced  by  the 
appropriate  term  from  the  equations  below. 


*  0 

e3  -  -  3/2(1  - e2)sin  ic oso>(  1  -  5/ 4  sin2!) 
e,  =  -45/16(1  -  e2)(  1  -7/6sin2i)esin2isin2ou 


=  15/4(1  - e2)sin (  (1  -  7  /2sin2 i*  21  /8sin 4 1)  ( 1  ■* 


7/8(1  9/8sin2Oe2sin2icos3u> 


] 


Pj  =  525 / 32 ( 1  -  e  )sin2i 


( 1  -  3  sin 2 1*  33/16  sin4  i)(  1 


3/16(1  1  1  /  1 0  sin  2 1  )e3  sin  2 1  sin 


3/4ez)cosuo 


+  1  /2ez)esin  2u> 


h  -  0 

i -  3/2(1  -  5/4sin20®cosoocosi 
1 4  -  45/32(1  -  7/6sin2i)e2sin2uusin2( 


u  =■  ■■  1  5/4erosi 


(1  -  7/2  sin2! +21/8  sin4i)(]  *  3 /4c2) cos co 
7/8(1  -  9/8sin2t)e2sin2icos3uu  j 
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£ 6  “  ~  525 /64c  sin  2 i (  1  -  3sin2 i  *  33/  1 6sin4 i) ( 1  -  1  /2e2)esin  2  to 
*3/16(1-11/  10sin2£)e3sin2£sin4toj 

-  ~  3/2cosi 

n3  =  3/2(1  -  15/4sin2£)esma>coti 

04  -  15/4 cos £( (  1  -7/4 sin2 0(1  +  3/2e2)-  3/4(1  -  7/3  sin  2  Oe*  cos  2  u>] 
fir  =  -  15/ 4 cot  i  (1-21  / 2  sin2  i  ■*  1 05  /8sin4  0(  1  *  3  /  4e2)c  sin  to 

*7/8(1  1  5/8sin2i)e3sin  3to 

Df  —  105/  16  cos  Jo  -  9/2  sin 2  i  +  33/8sin4 1)0  *  5e2*  15/8e4) 

■  5/2(1  6 sin 2 £  *  99/16 sin 4 0 (  1  *  1  /2e2)e2cos2io 

-  15/32(1  -  33/20sin2i)e4sin2£Cos4to 

a>t,  =  3(1  5  / 4  sin 2 1) 

to,  -  3/  2e  1  sin  tosin  £[  (  1  -5/4sinzi)*  (35 /4  cos2 1  -  cosec2 i)e2] 

u>,  -  1  5/32^  (  1  6  -  62sin2 1  *  49  sm4  £) -*  (6  sin2i  -  7  sin2  £)cos2u> ^  1  8 

-  63sin2i*  189/4sin4  £  jo2  *  (-6  *  35sin2£  -  63/2sin4£)e2cos2u> 
cor,--  1  05/ 1 6e' 1  sin  tocosec  i  ( -  4/7  *  2  sin2i  -  3/2sin 4  i)sin2  £  +  ^4/7 

-  87/7sin2£*  67/2 sin4i  -  357 /  I6sin‘i  je2  +  (  -  1  +  9/8 sin2 1) 

e 2  sin4  £  cos  2  to  *  (3/7-7sin2£  +  267  /  1 6sin4  £  -  165/16  sin6i)e4 
*  ( 1  -  39/8sinz£*  33/8sin40c4sin2*c°s2a>] 
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uo6  =  525/64  8/5(1  -  8sin2 1 1 29/8  sin4 1  -  297 /32  sin6  i)  -  ^2-6sin2i 

-*  33/8sm4 i  j sin  2 z cos 2uu  -*  6^1  -  43/6sin2  i  -*  109/8 si n4  z  -  121/8 

sin6  z  je2  +  (-2-*-25sin2z-  459/8  sin4  i  *  56 1  /  16sin6z)e2cos2u> 

•*3/8(1-11/  10sin2i)G2sin4 icos4au  +  ^2-27/2 sin2!-*-  99 / 4 sin 4 z 

-  429/32  sm  6  z)e4-*  (-  1  +  21  /2sin2i  -  363/  16sin4z  +  429/32 sin6 z) 

p4cog2ou  •*  3/8(-  1  22 /S  sin  2  z  -  143/40sin4z)e4sin2zcos4u> 

These  equations  may  be  used  to  determine  the  secular  changes  to  the 
orbited  elements  due  to  geopotential  which  occur  in  one  orbit  —  periodic 
perturbations  with  a  period  of  less  than  one  orbital  period  are  averaged 

out. 
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Air  Drag.  The  motion  of  a  satellite  through  the  upper  atmosphere 
of  a  planet  subjects  the  satellite  to  aerodynamic  forces  which  can 
perturb  the  orbit.  Atmospheric  forces  cause  an  acceleration  in  the 
direction  opposite  the  velocity  with  respect  to  the  atmosphere  (commonly 
called  drag),  and  an  acceleration  perpendicular  to  the  velocity  vector. 
The  perpendicular  component  does  not  necessarily  pass  through  the 
satellite's  center  of  mass,  resulting  in  both  a  "lift"  force  and  a  "lift" 
moment  about  the  center  of  mass.  Figure  1  illustrates  the  aerodynamic 
forces  on  a  satellite. 


Figure  1.  Atmospheric  Forces  on  the  Satellite 


The  direction  of  the  perpendicular  component  at  any  time  is  a  function 
of  satellite  attitude.  It  is  reasonable  to  assume  that  the  perpendicular 
component  is  small  compared  to  the  other  component,  and  that  the 
design  and  attitude  of  the  satellite  over  time  will  drive  the  resultant 
acceleration  towards  zero,  so  that  this  component  can  be  ignored  [8:  13]. 
The  acceleration  due  to  drag  is  then  given  by: 


CdA 

=  — p,v 


(16) 


where  C D  is  the  coefficient  of  drag  for  the  satellite,  A  is  the  frontal 
area  the  satellite  presents  to  the  air  stream,  m  is  the  satellite  mass,  p  is 
the  atmospheric  density,  and  v  is  the  velocity  with  respect  to  the  atmo¬ 
sphere.  To  simplify  parametric  analysis,  a  "ballistic  coefficient"  ft  is 
often  introduced  where  (3  -  Each  of  these  factors  will  now  be 
examined  in  turn. 

Ballistic  Coefficient.  Starting  with  ihe  simplest  term  of  the  ballis¬ 
tic  coefficient,  the  mass  of  the  satellite  may  be  either  fixed  or  a  func¬ 
tion  of  time.  The  mass  obviously  has  an  affect  on  the  acceleration  that 
drag  will  produce  on  the  satellite  (a  less  massive  satellite  will  be  more 
susceptible  to  drag  perturbations  than  a  more  massive  satellite)  in 
accordance  with  Newton's  Second  Law.  An  effective  or  average  value 
for  the  mass  over  the  period  of  interest  can  be  used  as  long  as  the 
actual  mass  does  not  vary  greatly  from  this  value. 

The  area  the  satellite  presents  to  the  atmosphere  depends  on  the 
size,  shape  and  orientation  of  the  satellite.  Uncontrolled  tumbling  or 
controlled  attitude  changes  will  cause  a  change  in  the  presented  area. 
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either  of  these  can  be  complex  functions  of  time  and/or  orbital  position. 
Again,  an  average  or  effective  value  must  be  employed. 

Due  to  the  long  mean  free  path  lengths  associated  with  the  molecu¬ 
lar  densities  encountered  at  orbital  altitudes,  Newtonian  flow  theory  is 
used  in  determining  the  coefficient  of  drag.  In  this  model  of  gas 
dynamics  the  shock  waves  are  considered  to  remain  very  close  to  the 
leading  edge  of  the  satellite,  so  each  molecule  acts  independently,  with 
nc  "knowledge"  of  any  other  molecule.  Several  major  assumptions  in 
this  model  are: 

1)  The  satellite  is  assumed  to  be  at  rest,  with  the  molecules  flow¬ 
ing  past  with  some  relative  velocity  plus  some  distribution  of  a  thermal 
velocity. 

2)  Each  molecule  will  impinge  upon  the  satellite  and  be  momentarily 
retained  before  being  re-emitted. 

3)  Collisions  between  incident  and  re-emitted  molecules  are 
neglected.  Some  broad  assumptions  must  be  made  in  evaluating  this 
coefficient,  but  several  sources  recommend  setting  2.0  <  C  D<  2.2  [8:  15;  9: 
2-4], 

Atmospheric  Density.  The  single  most  important  characteristic  of 
the  atmosphere's  role  in  determining  satellite  lifetime  is  density.  The 
atmospheric  density  is  a  function  of  time,  solar  activity,  and  altitude. 

As  the  planet  revolves  about  its  axis  the  sun  warms  the  atmosphere, 
causing  it  to  expand  and  form  a  "diurnal  bulge"  about  two  hours  after 
midday.  Although  this  causes  relatively  little  change  in  density  at  low 
altitudes,  density  can  increase  by  a  factor  of  8  at  higher  altitudes.  The 
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27  day  rotation  period  of  the  sun  causes  another  periodic  effect  on  the 
density.  The  position  of  the  planet  on  its  orbit  provides  another  oscil¬ 
lation,  with  minima  at  apohelion  and  maxima  at  perihelion.  Solar  activity 
is  related  to  the  11  year  sunspot  cycle,  and  this  also  affects 
atmospheric  heating.  Solar  activity  also  varies  essentially  randomly  from 
day  to  day  with  short  term  flares  and  other  activity. 

The  primary  assumption  made  in  modelling  the  atmosphere,  and  the 
greatest  simplification,  is  that  density  is  a  function  only  of  altitude 
above  some  reference  ellipsoid  (geodetic  altitude).  While  the  factors 
mentioned  in  the  preceding  paragraph  all  cause  deviations  from  this 
model,  empirical  data  from  actual  satellite  trajectories  indicates  density 
can  be  modelled  as  a  function  only  of  geodetic  altitude  reasonably  accu¬ 
rately  [14:  2;  8:  22].  Density  drops  off  with  altitude  in  accordance  with 
the  perfect  gas  law: 

P~4T  (17) 

(where  P  is  the  pressure,  M  is  the  mean  molecular  weight  of  the  atmo¬ 
sphere,  T  is  temperature  and  R'  is  the  universal  gas  constant)  and  the 
hydrostatic  equation: 

dp  -  -pgrdr  (18) 

Extensive  studies  of  earth  satellites  have  shown  a  good  first 
approximation  of  the  density  can  be  obtained  by  use  of  an  exponential 
atmosphere  model.  This  atmospheric  model  assumes  the  density  drops 
off  exponentially  with  altitude  h: 


where  p0  is  the  density  at  some  reference  altitude  h0,  and  H  is  the  scale 
height  (the  distance  over  which  density  will  change  by  a  factor  of  e). 
Several  types  of  exponential  atmospheres  exist,  with  the  primary  differ¬ 
ence  being  in  the  determination  of  the  scale  height  at  the  altitude  of 
interest.  Some  of  the  methods  of  determining  scale  height  are: 

1)  Strictly  Exponential  Atmosphere:  The  scale  height  is  be  consid¬ 
ered  to  be  constant  over  all  altitudes. 

2)  Locally  Exponential  Atmosphere:  The  scale  height  is  considered 
to  be  constant  over  some  small  altitude  interval. 

3)  (3r  Constant  Atmosphere:  The  product  of  the  inverse  of  the  scale 
height  (ft  not  to  be  confused  with  ballistic  coefficient)  and  the  radius 
considered  constant. 

4)  Isothermal  atmosphere:  The  temperature  is  considered  constant 
over  some  altitude  interval.  This  leads  to  the  quantity  £r2  being  con¬ 
stant  [14:  4], 

Each  of  these  methods  is  preferable  in  some  situations.  This  model 
uses  a  constant  scale  height  atmosphere,  with  the  scale  height  set  to 
the  proper  value  for  the  region  in  which  the  satellite  is  expected  to 
operate.  The  actual  scale  height  used  is  22.48  km  at  a  reference  alti¬ 
tude  of  250  km. 

Ellipsoidal  Altitude.  Since  the  density  is  a  function  of  ellipsoidal 
altitude,  it  is  necessary  to  find  ellipsoidal  altitude  as  a  function  of  lati¬ 
tude.  Given  a  planet  with  ellipticity  e,  the  ellipsoidal  altitude,  ha>  can 
be  approximated  by: 


where  <t>  is  the  geocentric  latitude  [9:  3-5]. 

Relative  Velocity.  The  velocity  v  the  satellite  has  with  respect  to 
the  rotating  atmosphere  is  the  difference  between  the  inertial  satellite 
velocity  (v,)and  the  velocity  of  the  atmosphere  at  that  point  (va): 

v=v.-va  (21) 


Finding  the  velocity  of  the  atmosphere  involves  several  assumptions. 

The  atmosphere  rotates  at  approximately  the  same  rate  as  the  planet 
near  the  surface,  but  this  rate  drops  off  with  altitude.  Assuming  that 
the  entire  atmosphere  rotates  uniformly  and  at  the  same  rate  as  the 
planet  greatly  simplifies  matters.  If  the  atmosphere  is  considered  to 
rotate  uniformly,  and  at  the  same  rate  as  the  planet,  then  the  inertial 
velocity  of  the  atmosphere  at  a  point  located  at  rais  given  by: 

v„  =  coa  x  r0  (22) 


with  uoa  being  the  inertial  rotation  rate  of  the  planet. 

Tins  model  includes  perturbations  due  to  atmospheric  drag.  In  fact, 
drag  is  the  probably  the  single  most  significant  factor  in  determining 
the  lifetime  of  low  orbiting  satellites,  and  orbits  for  which  drag  is  not 
the  predominant  perturbation  will  be  considered  beyond  the  scope  of 
this  work.  The  physical  parameters  of  the  satellite  used  in  this  model 
are  presented  in  Appendix  B. 

Averaged  Equations  for  Atmospheric  Drag.  Equations  representing 
the  average  change  in  orbital  elements  over  one  period  due  to  air  drag 
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have  been  produced  by  Sterne  [13].  These  are  presented  in  Equation 
(23),  where  E  is  the  eccentric  anomaly,  E  i  -  (  1  -  ecos£),  £z  -  (  1  -  e cos£),  n 
is  the  mean  motion  ((fi/a3)1/z),  and  d  -  (u>0/n)cosi(  1  -  e2)1'2. 


A«-2Po-Jppi-dL;) 


x  [cos  E  -  .Sd(£  ,  )(2cos  £-e-ecos2£)/(l  -  e2)]d£ 

2n 

A i  =•  -  .S(o)0/n)Pasini(l  -e2)'  S  J  p£*/2£j/2f  1  -  d^J  (23) 


x  <  1  -*  cos2uo£j2[(2  -  ez)cosz£  -  1  -*  2ez  -  2ecos£])d£ 


2* 

V 


AO  -  -  .!3(uja/n)Pasin2u)(]  -  e 2 )  "  I  p(  1  -  e2cos2£)l  1-d 


x  [2e2  -  1  -  2ecos£  +  (2  -  e2)cos2 £]d£ 


A oo  -  -  cos i AO 


These  equations  can  be  numerically  integrated  to  find  the  average 
change  in  the  orbital  elements  over  one  period  due  to  atmospheric  drag. 
As  with  the  geopotential  averaged  equations,  periodic  effects  will  be 
averaged  out. 
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Third  Body  Geometry.  If  a  third  body  is  close  enough,  or  massive 
enough,  its  gravitational  attraction  can  have  an  effect  on  the  satellite's 
orbit.  In  cases  where  the  third  body  is  in  a  nearly  circular  orbit  about 
the  primary,  or  the  primary  is  in  a  nearly  circular  orbit  around  the 
third  body,  the  third  body  can  be  modelled  as  being  in  two-body  motion 
around  the  planet,  simplifying  the  process  of  finding  the  vector  from 
the  third  body  to  the  planet. 


Figure  2.  Third  Body  Geometry 


Figure  2  shows  the  various  position  vectors  involved.  As  indicated, 
the  vector  from  the  planet  to  the  third  body  is  given  by  pIt>,  and  the 
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vector  r  locates  the  satellite  from  the  center  of  mass  of  the  planet. 
Hence,  the  vector  from  the  third  body  to  the  satellite  is  given  by 


P,  -  r  -  plt . 

Gravitational  Attraction.  The  effect  a  third  body  has  on  the  orbit 
of  a  satellite  is  two-fold:  1)  the  satellite  is  affected  by  the  attraction  of 
the  third  body,  perturbing  its  motion  around  the  planet,  and  2)  the 
motion  of  the  planet  itself  is  affected.  The  acceleration  due  to  the 
gravitational  attraction  of  the  third  body  in  both  cases  is  given  by: 


(24) 


where  the  vector  from  the  third  body  to  the  body  of  interest  is  given 
by  pie  and  n  lt  is  the  product  of  the  mass  of  the  third  body  and  the 
universal  gravitational  parameter.  The  complete  third  body  acceleration 
(including  the  direct  effect  of  the  third  body  on  the  satellite  and  the 
effect  on  the  planet)  is: 


f  P, 

Prt  ^ 

( (r , 

p^ 

a»  “  P.t, 

Ipf 

pH , 

- 

Ur“p»)'3 

pH  ) 

(25) 


This  model  includes  third  body  effects  from  the  sun  [9:  3-4], 

Averaged  Equations  for  Third  Body  Effects.  Kaufman  [6]  has 
developed  equations  representing  the  change  in  orbital  elements  over 
one  period  due  to  a  third  body.  In  the  event  that  the  inclination  is 
equal  to  zero  the  equations  for  the  change  in  inclination  and  longitude 
of  ascending  node  are  not  valid.  The  subscript  (©)  indicates  a  parame¬ 
ter  related  to  the  third  body.  These  equations  are  presented  below. 
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An  =  0 


A e  =  (-  lSn)(  1  -e2)  5(n./n)2eaP(a./r.)3 


Ai--6n(n0/n)2(l  -  e2)"  S(sint)  1  (a0/r9)3 
x[a6(l  +  4e2)  -*  (3c(l  -  e2)  -  5  ape 2  cost] 

AD-6n(nG/n)2(]-e?)'S(sint)'1(a0/r8)3  (26) 

x[a\(l  -  4e2)  sinai-*  py(  1  ~  e2)cosu>] 

Aa>  -  6n(n0/n)2(l  -  e2)"  S(a0/re)3x  |(4a2-p2-  1) 

-  (ycos  t  /  sin  t)(  1  e2)[a(i  -»  4e2)sinco  +  p(l  -  e2)cosco]l 


The  terms  a,  p,  v,  6  and  t  are  coordinate  frame  conversion  factors 

defined  by: 

a  -  r 0°  •  P 

p-  u°-q 

\  =  r0°-R 
6  -  r0° • P 

c  -  re°’Q' 
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As  before,  these  equations  can  be  used  to  obtain  changes  in  the  ele¬ 
ments  over  one  period  due  to  the  third  body. 
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Solar  Radiation  Pressure. 


Geometry.  Any  time  the  satellite  is  illuminated  by  the  sun,  solar 
radiation  pressure  exerts  a  slight  but  continuous  force  on  it,  producing 
a  small  acceleration.  This  is  caused  by  the  continuous  momentum 
transfer  from  innumerable  photon  impacts.  The  acceleration  depends  on 
the  area  of  the  satellite,  the  mass  of  the  satellite,  the  reflective 
characteristics  of  the  satellite  and  the  radiation  pressure  exerted  by  the 
sun  at  that  distance.  Figure  3  illustrates  the  various  vector  relations 
involved. 


Figure  3.  Solar  Radiation  Pressure  Geometry 


If  the  vector  from  the  planet  to  the  sun  is  ppG  and  the  vector  from  the 
planet  to  the  satellite  is  r,  then  the  vector  from  the  sun  to  the  satellite 
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is  Po,  *  r-pp0. 

The  radiation  pressure  exerted  is  a  function  of  distance  from  the 
sun  alone,  and,  like  gravity,  is  ruled  by  an  inverse  square  law.  Since 
the  radiation  pressure  on  a  perfectly  absorbing  (black  body)  surface  at 
the  distance  of  the  earth's  orbit  is  known,  it  is  convenient  to  express 
the  absorptivity  of  the  satellite  as  a  ratio,  and  the  distance  from  the 
satellite  to  the  sun  as  a  fraction  of  the  distance  from  the  earth  to  the 
sun. 

The  mean  distance  from  the  earth  to  the  sun  is  1  All,  or  149599650 
km.  Therefore,  the  solar  radiation  pressure  at  any  other  distance  from 
the  sun  will  be  related  to  the  pressure  at  1  AU  by  the  square  of  the 
product  of  the  inverse  of  the  distance  and  1  AU.  If  the  radiation  pres¬ 
sure  en  a  black  body  at  1  AU  is  Ps,  then  the  pressure  exerted  on  a 
black  body  at  the  satellite's  position  is 

Reflectivity  and  Other  Considerations.  In  reality,  no  satellite  (or 
any  other  object)  is  a  perfect  black  body.  The  reflectivity  of  the  satel¬ 
lite  car.  profoundly  change  the  magnitude  of  the  acceleration  due  to 
radiation  pressure.  The  acceleration  on  a  perfect  reflector  will  be  twice 
that  on  a  perfect  absorber,  and  a  transparent  satellite  would  feel  no 
acceleration  at  all.  Therefore  a  factor  v  is  included  in  the  expression 
for  the  acceleration,  where  v  “  1  for  a  perfect  absorber,  y  “  2  for  a  per¬ 
fect  reflector,  and  y  <  1  for  a  translucent  satellite,  with  read  world  satel¬ 
lites  usually  represented  by  values  between  1  and  2.  As  an  example, 
for  aluminum  y  -  1.96.  The  angle  of  reflection  can  adfect  the  direction  in 
which  the  acceleration  is  acting,  but  for  simplicity  it  is  justifiable  to 
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assume  the  satellite  reflects  diffusely,  resulting  in  the  force  pointing 
directly  away  from  the  sun. 

The  complete  description  of  the  acceleration  due  to  solar  radiation 
pressure  is  now: 


a.rp  “  vP/ 


A  AU 


m\ |pQ,| 


P*. 


(27) 


where  A  is  the  area  illuminated,  and  m  is  again  the  satellite  mass  [1: 


188]. 


Figure  4.  Occultation  Geometry 


Occultatron.  Note  this  acceleration  only  acts  when  the  sun  is  illu¬ 
minating  the  satellite  —  during  periods  when  the  satellite  is  eclipsed  or 
occulted  it  goes  to  zero.  To  check  for  planetary  occultation,  a  planet 
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centered  coordinate  system  is  constructed  with  the  first  unit  vector  (sj 
pointing  in  the  direction  of  the  sun  (Figure  4).  If  the  shadow  cast  by 
the  planet  is  assumed  to  be  a  cylinder,  and  the  satellite  is  located  from 
the  planet  by  the  vector  r,  the  satellite  will  be  occulted  at  any  time  that 
the  following  two  conditions  are  both  met: 

s  ,  •  r  <  1 

(F-£2)2+  (F-s3)2^rLm  (28) 

where  r  atrr.  is  the  radius  of  the  planet  plus  the  altitude  to  which  the 
atmosphere  can  be  considered  to  block  solar  radiation.  Then  at  any 
time,  the  geometry  of  the  sun-planet-satellite  system  can  be  checked  for 
the  presence  of  radiation  pressure  [10:  3-33],  Solar  radiation  pressure 
will  be  included  in  this  model. 

Averaged  Equations  for  Solar  Radiation  Pressure.  Cook  [4]  has 
derived  equations  for  the  rates  of  change  of  the  orbital  elements  over 
one  period.  In  these  equations  the  effects  of  the  primary’s  shadow  are 

neglected. 


A  ci  -  0 


,  3(l-os),/2^ 

- Zna  ^  p 


A  i  - 


3U'ecosuo 
2na( 1  - e 2 ) 1 ' 2 


An  - 


314/eSin  uo 


2na{]  e2)12  sin  1 


(29) 
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Aoo  - - -  S p-  An 

2nae  p 

The  terms  S p,  Sp  and  W  are  components  of  the  solar  radiation  force,  and 


are  defined  by: 


S  p  -  a  ,rp^  cos2  ^  cos  (oo  +  Cl  -  A)*  sin2^cos(co  +  Cl*  A)  cos2^ 

2  t,  ,  t  1  ,  I 

+  cos  -cos(oo  -  Cl  *  A)  sin  -cos(oo-n-A)  sin  - 


[cos(oo-  A)  -  cos(oo  +  A)]sin  isinf; 


Tp^-a.%rp[  cos2^sin(oo  +  Cl  -  A)  +  sin2|sin(co  +  Cl  *  A)  cos2^ 


cos  ^ sin  (oo  -  Cl  ->  A)  +  sin  -sin 


J  2 


(co-fl-A)  sin2^ 


[sin(oo-A)-sin(uo  +  A  )]sin  tsin  ^ 


t’  sin  oo  =: 


Cl  *  r  p  ! 


[cos (to  +  Cl-  A)  -  cos  (uo  ~  n+  A)]  sin  ( cos' 


[cos(uo  +  H  +  A  )  -  cos(oo  -  f)-A)]sinisin‘ 


-  [cos(ou  *  A)-cos(oo  -  A  )]cosi  sin  t. 


cos  oo  —  y[sm(co  -  n-A)-sin(oo-D  +  A)]sini  cos2  - 


[sin(to  +  n  +  A)  -  sin(ou  -  Cl  -  A)]sinisin  - 


[sin(oo  A)-sin(oo  -  A)]cosisin£ 


where  a,rp  is  the  magnitude  of  the  acceleration  due  to  solar  radiation 
pressure  described  by  Equation  (27),  t,  is  the  obliquity  of  the  ecliptic, 
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and  A  is  the  geometric  mean  longitude  of  the  sun,  measured  in  the 
ecliptic.  These  equations  may  be  used  to  find  the  secular  and  long 
period  perturbations  to  the  orbit  due  to  solar  radiation  pressure. 
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Numerical  Model 


Two  methods  of  applying  all  of  the  aforementioned  perturbations  in 
a  mathematical  model  of  satellite  motion  were  considered  —  an  averaging 
method  and  Cowell's  method.  An  averaging  method  would  use  the 
averaged  equations  presented  at  the  end  of  each  of  the  preceding 
sections  (Equations  (15),  (23),  (26)  and  (29))  to  propagate  the  average 
orbital  elements  foruard  in  time,  giving  a  prediction  of  what  the 
elements  would  be  at  the  end  of  the  90  day  period.  Cowell's  method 
would  propagate  the  position  and  velocity  of  the  satellite  forward  in 
time,  subject  to  equations  of  motion  found  by  summing  the  acceleration 
equations  presented  in  each  section  (Equations  (13),  (16),  (25)  and  (27)). 
Both  methods  have  advantages,  and  both  were  investigated  in  making 
the  selection. 

Cowell's  method  of  General  Perturbations  was  chosen.  In  this 
method,  the  equations  of  motion  for  the  satellite,  including  the 
perturbative  accelerations,  are  expressed  in  rectangular  coordinates  and 
integrated  numerically  [12:  199].  This  allows  the  use  of  the  acceleration 
equations  determined  in  each  section,  without  requiring  the  averaged 
equations.  The  starting  orbital  elements  are  converted  to  a  state  vector 
in  rectangular  coordinates  and  used  as  initial  conditions  for  the 
integration  of  the  equations  of  motion.  At  each  time  step,  the  state 
vector  is  converted  back  to  orbital  elements  so  that  they  may  be 
examined,  if  desired.  This  method  has  the  advantages  of  being  simple 
and  straightforward,  in  addition  to  determining  the  actual  position  of 
the  satellite  at  each  time  step.  Its  disadvantage  is  that  it  is  slower 
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than  methods  which  work  with  averaged  orbital  elements. 
The  equations  for  this  model  are  of  the  form: 


x  = 

y  =  t\ 

z  =  vy 


X  -  vx 

y  ~v> 

z  -  u. 


^(J1,  a<*  X  *  +  a  irp  x 

a  fy*  a  “y*  a,t  y  *  a  *rPy 
a  a  z  *  a  «  i  "*  ^tbjr  +  a  trp  2 


(30) 


where  a a.,  u..,  and  as.p  are  given  in  Equations  (13),  (16),  (25),  and 
[Z~"\  respectively.  The  motion  is  expressed  in  six  first  order  equations 
to  allow  use  of  a  fixed  step  size,  eighth  order  Runge-Kutta  integration 
algorithm  described  in  reference  [5]. 


m.  Results 


Search  Method 

Orbital  elements  on  the  Critical  Surface  were  obtained  by  estimating 
an  initial  element  set,  propagating  that  orbit  forward  90  days,  and  itera¬ 
tively  adjusting  the  estimate  and  repropagating  the  orbit  to  obtain  a 
final  value  of  periapsis  altitude  within  a  given  tolerance  of  the  decay 
threshold.  This  method  was  necessary  because  the  numerical  integration 
routine  used  was  incapable  of  being  operated  backwards  in  time. 

It  was  first  necessary  to  determine  realistic  decay  criteria.  A 
decay  threshold  chosen  too  low  resulted  in  very  rapid  decay  at  the  90 
day  point,  making  it  very  difficult  to  "fine  tune"  the  initial  semi-major 
axis  such  that  on  day  90  the  periapsis  altitude  would  be  within  a  given 
tolerance  of  the  decay  threshold.  If  the  decay  threshold  was  too  high, 
drag  would  not  be  the  predominant  perturbation,  and  the  orbit  would 
tend  to  not  decay  monotonically,  with  perturbations  other  than  drag 
having  a  more  pronounced  affect  on  satellite  lifetime. 

Once  a  suitable  threshold  was  determined,  the  initial  estimate  was 
obtained  by  computing  the  two-body  semi-major  axis  for  the  desired 
periapsis  altitude.  Recalling  Equation  (3)  from  Chapter  II: 


and  recognizing  that  at  perigee  v  -  Q,  the  radius  at  perigee  can  be 
expressed  as: 

rp-a(l-e)  (32) 

This  initial  estimate  was  then  propagated  forward  through  90  days 
using  the  numerical  method  described  in  Chapter  II.  If  the  final  per- 
iapsis  altitude  was  not  within  the  decay  threshold  tolerance,  the  semi¬ 
major  axis  was  incrementally  adjusted  until  the  desired  periapsis 
altitude  was  obtained.  This  procedure  was  repeated  over  a  range  of 
inclinations  at  each  of  five  eccentricities.  Due  to  the  iterative  nature  of 
the  search  method,  obtaining  a  single  data  point  required  between  ten 
and  fifteen,  and  sometimes  as  many  as  twenty,  integrator  runs.  Since 
each  90  day  integrator  run  required  approximately  twenty  minutes  of 
computer  time,  producing  the  Critical  Curves  was  a  very  time  intensive 
operation. 

A  polynomial  fitting  routine  was  then  applied  to  the  data  to  gener¬ 
ate  functions  which  would  allow  interpolation  between  the  data  points 
and  simplify  use  of  the  Critical  Curve  at  that  eccentricity. 
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D&mi  Threshold 

To  establish  a  realistic  decay  threshold,  general  decay  characteris¬ 
tics  were  examined  at  e  -  .001.  For  a  given  periapsis  altitude,  a  less 
eccentric  orbit  will  have  a  smaller  semi- major  axis  than  a  more  eccentric 
orbit,  and  will  spend  a  larger  fraction  of  the  orbital  period  at  lower 
altitudes  where  the  satellite  will  be  subjected  to  drag.  Therefore,  lower 
eccentricity  orbits  are  more  profoundly  affected  by  drag,  so  a  realistic 
threshold  for  low  eccentricity  should  be  acceptable  for  higher  eccentri¬ 
cities.  The  behavior  of  the  periapsis  altitude  over  time  for  a  typical 
orbit  (a,  -  6270km ,  e,  -  .001  ,  i,  -  45  )  is  shown  in  Figure  5. 


I 


Figure  5.  Terminal  Orbital  Decay 
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It  can  be  seen  that  after  a  periapsis  altitude  of  150  km  is  reached, 
decay  becomes  very  rapid.  At  periapsis  altitudes  below  about  120  km, 
the  curve  is  almost  vertical.  Attempting  to  determine  initial  conditions 
to  result  in  a  final  periapsis  altitude  in  this  region  would  not  only  be 
very  difficult,  but  the  lifetime  of  any  orbit  in  that  regime  would  be 
extremely  sensitive  to  variations  in  initial  conditions  and  perturbative 
forces.  Therefore,  a  minimum  decay  threshold  of  130  km  was  selected. 
Two  additional  decay  threshold  altitudes  were  selected  at  140  km  and 
150  km.  The  tolerance  allowed  was  10  km  above  the  threshold.  When 
the  slope  of  the  decay  curve  permitted,  an  effort  would  be  made  to 
obtain  a  final  periapsis  altitude  as  close  to  the  threshold  as  possible. 
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Critical  Cnrv» 

Eccentricity  =  .001  The  smallest  eccentricity  examined  was  e-  .OOl 
to  reflect  the  scarcity  of  perfectly  circular  orbits  in  real  world 
situations.  The  curve  family  is  shown  in  Figure  6. 


Figure  6.  Critical  Curve  for  e=.001 

A  high  degree  of  inclination  dependence  can  be  seen  in  these  curves, 
causing  a  variation  in  initial  semi-major  axis  of  nearly  5  km  over  the 
inclination  range.  This  seems  to  indicate  a  much  lower  rate  of  decay  at 
low  and  high  inclinations  than  at  mid-inclinations.  A  short  investigation 
as  to  the  nature  and  cause  of  this  inclination  dependence  is  included 
later  in  this  chapter. 
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The  polynomial  fitting  routine  was  unable  to  produce  a  good  fit  to 
the  data  with  less  than  a  seventh  order  polynomial,  reflecting  the  rela¬ 
tive  complexity  of  the  curve.  The  polynomial  coefficients  are  given  in 


Appendix  C,  Table  C-l. 

Eccentricity  =  .01  At  e-  .01  the  inclination  dependence  is  still 
obvious,  although  less  pronounced.  The  curve  family  is  shown  in  Figure 
7. 


Figure  7.  Critical  Curve  for  e=.01 


Although  the  inclination  dependence  produces  a  less  striking  effect  in 
the  shape  of  the  curve  at  this  eccentricity,  the  effect  on  the  range  of 
initial  semi- major  axis  values  is  very  noticeable,  with  almost  25  km  of 
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variation . 


The  polynomial  fitting  routine  produced  a  good  fit  to  these  data 
sets,  and  all  of  the  following  data  sets,  with  a  fourth  order  polynomial. 
The  polynomial  coefficients  are  given  in  Appendix  C,  Table  C-2. 

Eccentricity  =  .02  At  e  -  .02  it  was  first  seen  that  at  higher  incli¬ 
nations  (as  before,  corresponding  to  slower  decay)  periapsis  altitude 
actually  increased  over  the  90  day  period. 


Figure  8.  Critical  Curve  for  e=.02 

Data  showing  this  tendency  is  not  presented  here,  since  it  is  not  of 
interest  in  terms  of  short  term  drag-limited  lifetime.  This  results  in 
some  of  the  fitted  curves  being  clipped  at  higher  inclinations.  The 
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curves  presented  here  are  clipped  at  the  first  inclination  where  the 
final  periapsis  altitude  is  higher  than  the  initial  periapsis  altitude.  The 
curve  family  is  shown  in  Figure  8.  The  trend  towards  increasing  range 
of  initial  semi-major  axis  values  continued,  spanning  almost  45  km  in 
this  case. 

The  polynomial  fit  was  again  accomplished  with  a  fourth  order  poly¬ 
nomial.  The  polynomial  coefficients  are  given  in  Appendix  C,  Table  C-3. 

Eccentricity  =  .04  At  e  -  .04  the  tendency  mentioned  in  the  last 
section  became  more  pronounced,  requiring  further  clipping  of  the 
curves.  The  curve  family  is  shown  in  Figure  9. 


Figure  9.  Critical  Curve  for  e=.04 


-39- 


The  trend  towards  a  wider  spread  in  initial  semi-major  axis  seemed  to 
have  leveled  off  at  just  over  50  km  in  this  and  the  next  graph,  but  this 
can  largely  be  attributed  to  the  data  clipping.  The  tendency  for  per- 
iapsis  altitude  to  "climb’'  at  high  inclinations  would  significantly  increase 
the  range,  were  the  higher  inclinations  to  be  included. 

The  polynomial  coefficients  for  these  curves  are  given  in  Appendix 
C,  Table  C-4. 

Eccentricity  =  .1  The  largest  eccentricity  examined  was  e=  .1.  The 
curve  family  is  shown  in  Figure  10. 


Figure  10.  Critical  Curve  for  e=.l 


It  is  obvious  that  at  this  and  larger  eccentricities,  smaller  inclina¬ 
tion  ranges  will  "decay  down”  to  a  fined  periapsis  adtitude.  At  these 
eccentricities,  drag  will  be  less  of  a  factor  in  determining  lifetime,  so 
these  will  no  longer  be  short  term,  low  orbiting  satellites. 

Again,  a  fourth  order  polynomial  was  sufficient.  The  polynomial 
coefficients  are  given  in  Appendix  C,  Table  C-5. 

Note  that  in  this  section,  the  x-axis  represents  initial  semi- major 
axis,  so  "peaks"  on  these  graphs  are  indicative  of  higher  amounts  of 
decay. 
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To  investigate  the  cause  of  the  inclination  dependence  noted  in  the 
graphs  above,  the  model  was  run  for  90  days  in  five  different 
configurations.  In  four  of  the  trials,  all  of  the  perturbations  but  one 
were  zeroed  out.  In  the  fifth  all  of  the  perturbations  were  zeroed  out. 
The  initial  conditions  used  were  considered  to  be  representative  of  the 
range  of  initial  conditions  considered  in  determining  the  Critical  Curves. 
The  five  curves  presented  in  these  graphs  represent  the  final  perigee 
altitude,  as  a  function  of  inclination,  for  each  of  the  separate 
perturbations  and  the  unperturbed  case.  The  no  perturbation  case 
differs  from  the  initial  perigee  altitude  by  only  the  integrator  error. 
These  graphs  depict  final  periapsis  altitudes,  so  the  "peaks"  represent 
minimal  decay,  the  opposite  of  the  previous  section. 

As  can  be  seen  from  Figure  11,  the  inclination  dependence  at  middle 
eccentricities  resides  almost  entirely  in  the  geopotential  perturbation. 

The  drag  perturbation  causes  the  same  total  decay  at  all  inclinations. 
This  is  different  than  would  be  the  case  for  the  earth,  where  the 
equatorial  bulge  in  the  planet  is  mirrored  in  an  equatorial  bulge  in  the 
atmosphere  (and  in  the  ellipsoidal  altitude  density  model  used  here)  and 
drag  perturbations  are  therefore  inclination  dependent,  with  low 
inclinations  being  more  affected  by  drag.  Because  Venus  is  much  more 
spherical  than  earth,  there  is  no  significant  corresponding  equatorial 
bulg®.  The  solar  third  body  perturbation  seems  to  be  negligible,  having 
less  of  an  effect  than  even  the  solar  radiation  pressure.  The  radiation 
pressure  perturbation  causes  a  slight  inclination  dependence,  with 


maximum  effect  near  20  degrees  of  inclination.  Neither  radiation 
pressure  nor  third  body  effects  are  normally  considered  to  be 
significant  for  low  orbiting  satellites,  but  were  included  here  because  of 
Venus'  proximity  to  the  sun. 

Figure  12  presents  the  same  data  as  Figure  11,  but  for  the  low 
eccentricity  case.  In  this  case,  the  same  observations  can  be  made 
concerning  solar  radiation  pressure  and  third  body  effects.  These  could 
be  expected  to  have  even  less  of  an  effect  in  this  lower  orbit.  The 
drag,  as  expected,  is  much  more  predominant  for  this  lower  orbit,  and  is 
again  independent  of  inclination.  The  geopotential  effect  is  significantly 
different,  with  almost  no  tendency  to  cause  the  orbit  to  "climb." 
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Figure  11.  Final  Periapsis  Altitude  as  a  Function  of  Inclination  (e=.02) 
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Figure  12.  Final  Periapsis  Altitude  as  a  Function  of  Inclination  (e-.OOl) 
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Since  the  geopotential  produces  a  conservative  force,  it  could  be 
expected  that  the  dependence  seen  in  Figures  11  and  12  would  be  part 
of  a  cyclic  behavior.  To  investigate  this  hypothesis,  longer  (120  day) 
simulations  at  the  previously  investigated  eccentricities  of  .02  and  .001 
were  run. 

High  Eccentricity.  Inclinations  used  for  the  higher  eccentricity 
case  were  40  degrees  (where  the  gravity  perturbation  leads  to  the 
fastest  decay),  65  degrees  (where  the  gravity  perturbation  has  a  minimal 
effect),  and  85  degrees  (where  the  gravity  perturbation  seems  to 
actually  make  the  periapsis  altitude  increase). 


Figure  13.  Periapsis  Altitude  as  a  Function  of  Time 
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Figure  13  shows  the  effect  geopotential  has  on  periapsis  altitude. 
As  could  be  expected  from  Figure  11,  this  shows  up  as  a  decrease  in 
periapsis  altitude.  The  expected  periodicity  (due  to  geopotential  pro¬ 
ducing  a  conservative  force)  is  not  apparent,  however.  Note  that  in 
these  graphs,  like  Figures  11  and  12,  the  final  periapsis  altitude  is 
plotted,  so  lower  parts  of  the  curve  correspond  to  more  decay. 


little  decay  at  65  degrees  inclination.  This  short  period  oscillation 
would  also  show  up  in  Figure  13  but  is  lost  on  that  graph  because  its 
magnitude  is  so  much  smaller  than  that  of  the  "net"  change. 


Figure  15.  Periapsis  Altitude  as  a  Function  of  Time 
In  Figure  15  the  expected  "climb"  in  periapsis  altitude  can  be  seen. 
All  three  of  these  figures  show  some  magnitude  of  net  change  over  the 
120  days,  along  with  some  small  fast  oscillation.  In  no  case  does  the 
larger  net  change  give  evidence  of  being  periodic. 
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Low  Eccentricity.  Inclinations  used  for  the  low  eccentricity  case 
were  40  degrees  (where  the  most  decay  occurred),  65  degrees  (where 
the  geopotential  had  minimal  effect),  and  70  degrees  (where  the  geopo¬ 
tential  caused  a  small  increase  in  periapsis  altitude).  The  results  are 
presented  in  Figures  16,  17,  and  18, 


Figure  16.  Periapsis  Altitude  as  a  Function  of  Time 
Figure  16  is  very  similar  to  Figure  13,  showing  the  overall  drop  in 
periapsis  altitude,  consistent  with  Figure  12. 
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Figure  17.  Periapsis  Altitude  as  a  Function  of  Time 
In  Figure  17  the  small  net  change  is  again  seen,  with  the  expected 
small  fast  oscillation  visible  because  of  the  scale  of  the  graph. 


;§.  2  1  7.00 

K 


5  2  1  6.00  ■ 


Figure  18.  Periapsis  Altitude  as  a  Function  of  Time 
In  Figure  18,  evidence  of  the  expected  periodicity  in  periapsis  alti¬ 
tude  is  finally  seen.  Because  this  curve  represents  the  relatively  small 
"climb''  in  periapsis  altitude  seen  at  low  eccentricity,  the  short  period 
oscillation  is  more  noticeable  than  in  Figure  15.  These  figures,  taken 
together,  show  that  the  effect  of  geopotential  seen  with  this  model  is 
periodic,  as  predicted  by  the  fact  that  the  geopotential  produces  a  con¬ 
servative  force.  The  response  is  composed  of  multiple  oscillations  of 
differing  periods.  This  investigation  implies  that  the  inclination  of  65 
degrees  is  relatively  stable  for  the  periapsis  altitude  for  Venus.  This  is 
accompanied  by  relative  stability  in  semi-major  axis  and  eccentricity, 
with  changes  in  these  elements  for  both  examined  eccentricities  of  less 
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than  .25  km  and  .0004,  respectively,  over  the  120  day  period.  Satellites 
placed  in  65  degree  inclination  orbits  would  experience  minimal  change 
in  periapsis  altitudes,  at  least  within  the  120  day  time  frame  examined. 
It  should  be  possible  to  adjust  the  inclination  so  that  geopotential 
induced  "climb"  will  offset  the  drag  induced  decay  in  the  periapsis,  at 
higher  altitudes,  resulting  in  short  term  orbits  with  roughly  constant 
periapsis  altitudes. 
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Visual  representations  of  the  three  Critical  Surfaces  were  produced 
with  the  aid  of  a  three-dimensional  plotting  program.  These  plots  are 
presented  as  Figures  19,  20  and  21. 

Note  that  in  these  plots,  inclination  runs  along  the  nearest  axis, 
increasing  from  left  to  right.  Eccentricity  (multiplied  by  1000  for 
scaling)  is  plotted  along  the  other  horizontal  axis,  increasing  from  front 
to  back.  The  vertical  axis  represents  the  difference  between  initial 
semi-major  axis  and  the  planetary  radius. 

The  Critical  Curves  can  be  seen  as  the  lines  of  constant 
eccentricity  corresponding  to  the  values  discussed  earlier.  Other  values 
on  the  surface  are  interpolations  of  the  plotting  program.  Obviously, 
more  data  would  significantly  increase  the  accuracy  of  these  surfaces, 
but  the  trend  of  the  data  is  apparent.  The  plotting  program  also 
interpolates  over  regions  where  the  data  has  been  clipped  at  higher 
eccentricities  and  inclinations. 

Any  set  of  orbital  elements  corresponding  to  a  point  above  these 
Critical  Surfaces  will  produce  an  orbit  which  will  not  decay  below  the 
specified  minimum  altitude  in  90  days. 


IV.  Conclusions  and  Recommendations 

Critical  Surfaces  can  be  produced  with  an  iterative  numerical 
method.  These  surfaces  show  a  strong  geopotential-caused  inclination 
dependence  which  is  periodic  in  nature.  Critical  Surfaces  generated  for 
lifetimes  in  excess  of  90  days  would  be  more  effected  by  this  long  term 
periodicity.  Drag,  as  expected,  causes  a  uniform  decay  rate  at  all  incli¬ 
nations  when  considered  for  a  spherical  planet.  Although  neither  solar 
third  body  effects  nor  radiation  pressure  made  substantial  contributions 
to  the  rates  of  decay,  solar  radiation  pressure  did  have  a  small  effect, 
and  would  probably  be  significant  for  higher  orbits  because  of  Venus' 
proximity  to  the  sun. 

At  certain  inclinations  and  eccentricities  the  periodic  geopotential 
effect  takes  the  form  of  an  increase  in  periapsis  altitude.  This  could  be 
used  to  offset  the  decay  due  to  drag  for  some  finite  time,  in  order  to 
maintain  a  constant  periapsis  altitude.  This  possibility  should  be  inves¬ 
tigated  more  closely.  For  high  eccentricity  and  inclination,  the  geopo- 
tential  effect  offset  the  drag  decay  completely,  resulting  in  orbits  which 
were  not  lifetime -limited  by  decay  in  the  time  frame  examined.  Since  the 
geopotential  effect  is  periodic,  these  orbits  would  have  finite  lifetimes, 
but  this  effect  could  be  used  in  conjunction  with  limited  station - 
keeping,  and  this  also  should  be  investigated. 

This  investigation  only  considered  zonal  harmonics  up  to  sixth 
order.  Further  investigation  involving  sectoral  and  tesseral  harmonics 
could  be  conducted,  but  would  require  consideration  of  all  six  orbital 
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elements,  complicating  analysis. 

As  mentioned  in  the  Introduction,  perhaps  the  most  interesting 
application  of  this  type  of  investigation  would  be  for  lightsat  mission 
planning.  Interplanetary  exploration  probes  are  not  likely  to  be  short¬ 
term  expendable  spacecraft  in  the  near  future,  but  the  military  potential 
of  small,  inexpensive,  single  purpose  satellites  is  great.  Producing  the 
Critical  Surfaces  for  these  types  of  missions  would  require  the  inclusion 
of  third  body  effects  from  the  moon,  which  would  add  complexity  to  the 
surfaces.  This  area  of  research  should  be  pursued. 


# 
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Appendix  A 


flmial  gomriaata  for  the  pi*™*  yhmi8 

The  physical  constants  for  the  primary  used  in  this  model  were 
taken  from  References  9  and  10. 


Mean  Distance  from  the  Sun 
Planetary  Radius  (i?,) 
Occultation  Radius  (ra(m) 

U. 

Rotation  Rate  (uoj 
Planetary  Ellipticity  (fc) 

Scale  Height  ((3) 

Reference  Height 
Reference  Density 


108100000  km 
6051  km 
6141  km 

3.24858  E  5  kmVsec2 
-1.71  E  -5  deg/sec 
0.0 

22.48  km 
250  km 

3.19  E  -4  kg/km3 


Geopotential  Coefficients 


J  2 

-.45207 

E  -5 

J  3 

.13421 

E  -5 

J * 

.24135 

E  -5 

J  S 

.25940 

E  -6 

.33613 

E  -6 
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Appendix  B 

Ehiaial  at  »V  spacecraft  tor  this  Model 

The  physical  parameters  for  the  spacecraft  for  this  model  were 
loosely  based  on  those  for  the  Magellan  Venus  Orbiter  (Reference  9). 


Spacecraft  Mass  (m) 

1085  kg 

Coefficient  of  Drag  (CJ 

2.0 

Spacecraft  Area  (A) 

24  m2 

Reflectivity  Factor  (y) 

1.0 
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The  coefficients  for  the  polynomials  fitted  to  the  Critical  Curve  data 
are  given  in  Table  C-l  through  C-5. 


Table  C-l.  Critical  Curve  Polynomial  Coefficients  for  e=.001 


DEGREE 


DECAY  THRESHOLD  (Periapsis  Altitude)  | 

130  km 

140  km 

150  km 

6267.99 

6268.24 

6268.49 

-6.95906  E-2 

2.99023  E-2 

1.32518  E-l 

2.08757  E-2 

7.38482  E-3 

-7.20412  E-3 

-1.28716  E-3 

-4.75775  E-4 

4.50068  E-4 

4.49841  E-5 

2.08839  E-5 

-6.96061  E-6 

-8.90591  E-7 

-5.18589  E-7 

-9.0554  E-8 

8.79975  E-9 

5.92692  E-9 

2.63286  E-9 

-3.33571  E-ll 

-2.45627  E-ll 

-1.44104  E-ll 

Table  C-2.  Critical  Curve  Polynomial  Coefficients  for  e=.01 


DEGREE 

DECAY  THRESHOLD  (Periapsis  Altitude)  | 

130  km 

140  km 

150  km 

0 

6287.83 

6288.54 

6289.48 

1 

1.77987  E-2 

-2.62089  E-2 

-9.69216  E-3 

2 

2.62497  E-2 

3.32471  E-2 

4.07476  E-2 

3 

-6.79421  E-4 

-8.40906  E-4 

-1.04491  E-3 

4 

4.20771  E-6 

5.20687  E-6 

6.53673  E-6 
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Table  C-3.  Critical  Curve  Polynomial  Coefficients  for  e=.02 


DEGREE 

DECAY  THRESHOLD  (Periapsis  Altitude)  j 

130  km 

140  km 

150  km 

0 

6331.34 

6333.85 

6337.94 

1 

2.75642  E-2 

-4.46667  E-3 

8.94402  E-2 

2 

4.31136  E-2 

5.70271  E-2 

6.14433  E-2 

3 

-1.11557  E-3 

-1.4776  E-3 

-1.61977  E-3 

4 

6.92373  E-6 

9.3038  E-6 

1.01088  E-5 

Table  C-4.  Critical  Curve  Polynomial  Coefficients  for  e=.04 


DEGREE 

DECAY  THRESHOLD  (Periapsis  Altitude)  | 

130  km 

140  km 

150  km 

0 

6447.23 

6454.59 

6463.24 

1 

3.03818  E-l 

4.0992  E-l 

5.80655  E-l 

2 

4.82529  E-2 

4.42407  E-2 

3.25586  E-2 

3 

-1.32737  E-3 

-1.23606  E-3 

-9.41725  E-4 

4 

7.95499  E-6 

7.12796  E-6 

4.68943  E-6 

Table  C- 

5.  Critical  Curve 

Polynomial  Coefficients  for  e=.l 

DEGREE 

DECAY  THRESHOLD  (Periapsis  Altitude)  f 

130  km 

140  km 

150  km 

0 

6875.97 

6882.02 

6891.97 

1 

-2.89251  E-l 

2.33786  E-l 

6.39414  E-l 

2 

7.58337  E-2 

5.74248  E-2 

2.77024  E-2 

3 

-1.87211  E-3 

-1.60948  E-3 

-9.46644  E-4 

4 

1.02871  E-5 

5.81199  E-6 
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